{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "if not any(path.endswith('textbook') for path in sys.path):\n",
    "    sys.path.append(os.path.abspath('../../..'))\n",
    "from textbook_utils import *\n",
    "\n",
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(sec:theory_randomAssignment)=\n",
    "# Simulating a Randomized Trial: Vaccine Efficacy\n",
    "\n",
    "In controlled experiments, we want to assign the participating individuals to one of two (or more) groups. In a drug trial, volunteers for the trial either receive the new treatment or a placebo (a fake treatment). In an A/B test of a new feature on a Web site, visitors to the site would either see the new feature or the usual Web page.  In both examples, we control the assignment to groups, and in a *randomized* controlled experiment, we use a chance process to make the assignment.  \n",
    "\n",
    "In drug trials, the scientists often use essentially an urn model model to select the subjects for the treatment, where those not selected receive the placebo. With A/B testing, we often use a systematic approach, where, for example, every other visitor to the page is shown the new feature (see the exercises to learn more about systematic sampling). \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We explore the SRS approach to random assignment with a small group of just $6$ individuals, labeled  $A - F.$ We enumerate all the possible treatment  groups that can be formed by a SRS of size $3$ from this population:\n",
    "\n",
    "$$ABC ~~ ABD ~~ ABE ~~ ABF \\\\ ACD ~~ACE ~~ ACF ~~ ADE \\\\ ADF ~~ AEF ~~ BCD ~~ BCE \\\\ BCF ~~ BDE ~~ BDF ~~ BEF \\\\ CDE ~~ CDF ~~ CEF ~~ DEF $$\n",
    "\n",
    "Each sample represents a possible treatment group, and its complement the control group, e.g., the first sample has $ABC$ in the treatment group, and $DEF$  in the control group. As we saw in {numref}`sec:theory_samplingVariation`, we can use the chance mechanism to work with probabilities of particular subsets and apply these to the outcome of an experiment. The next section provides an example.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Efficacy of COVID-19 Vaccines\n",
    "\n",
    "Detroit Mayor Mike Duggan made national news in March 2021 when he turned down a shipment of over 6,000 Johnson & Johnson vaccine doses stating that the citiziens of his city should \"get the best\". The mayor was refering to the efficacy rate of the vaccine, which was reported to be about 66%. In comparison, Moderna and Pfizer both reported efficacy rates of about 95% for their vaccines. \n",
    "\n",
    "Duggan's reasoning seems valid, but the scope of the three clinical trials are not comparable, meaning direct comparisons of the experiments is problematic {cite}`irfan2020`. Moreover, the Centers for Disease Control (CDC) considers a 66% efficicay rate quite good, which is why it was given emergency approval {cite}`cdc2021`.\n",
    "\n",
    "We consider each of these points in turn, beginning with scope and then efficacy.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scope\n",
    "\n",
    "Recall that when we evaluate the scope of the data, we consider who, when, and where the study took place. For the Johnson & Johnson clinical trial, participants: \n",
    "\n",
    "+ included adults 18 and over, where roughly 40% had conditions, called comorbidities, associated with an increased risk for getting severe COVID-19;\n",
    "+ enrolled in the study from October to November, 2020;\n",
    "+ were from 8 countries across 3 continents, including the US and South Africa.\n",
    "\n",
    "The participants in the Moderna and Pfizer trials were primarily from the US, roughly 40% had comorbidities for severe COVID-19, and the trial took place earlier, over summer 2021. \n",
    "The timing and location of the trials make them difficult to compare.\n",
    "Cases of COVID-19 were at a low point in the summer in the US, but they shot up in the late fall.\n",
    "Also, a variant of the virus that is more contagious was spreading rapidly in South Africa at the time of the J&J trial. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Each clinical trial was designed to test a vaccine against no vaccine under similar circumstances. \n",
    "This is accomplished through random assignment of subjects to treatment and control groups. \n",
    "While the scope from one trial to the next are quite different, \n",
    "the randomization within a trial keeps the scope of the treatment and control groups roughly the same, which enables meaningful comparisons between groups in the same trial. \n",
    "\n",
    "The scope was different enough across the three vaccine trials to make direct comparisons problematic. \n",
    "Next, we consider what the vaccine's efficacy rate means.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Design\n",
    "\n",
    "Johnson & Johnson enrolled 43,738 people in their trial {cite}`janssen2021`.\n",
    "These participants were split into two groups at random.\n",
    "Half received the new vaccine, and the other half received a placebo, such as a saline solution.\n",
    "Then, everyone was followed for 28 days to see whether they contracted COVID-19.\n",
    "\n",
    "A lot of information was recorded on each patient, such as their age, race, and sex, and in addition,\n",
    "whether they caught COVID, including the severity of the disease. At the end of 28 days, they found 468 cases of COVID-19, 117 of these were in the treatment group, and 351 in the control group.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "The random assignment to treatment and control, gives the scientists a framework to assess the effectiveness of the vaccine. \n",
    "The typical reasoning goes as follows:\n",
    "\n",
    "+ Begin with the assumption that the vaccine is ineffective\n",
    "+ so the 468 who caught COVID-19 would have caught it whether or not they received the vaccine\n",
    "+ and, the remaining 43,270 people in the trial who did not get sick would have remained healthy whether or not they received the vaccine.\n",
    "+ The split of 117 sick people in treatment and 351 in control was just do to the chance process in assigning participants to treatment or control.\n",
    "\n",
    "We can set up an urn model to reflect this scenario.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Urn Model\n",
    "\n",
    "Our urn has 43,738 marbles, one for each person in the clinical trial.\n",
    "Since there were 468 cases among them,\n",
    "we label 468 marbles with a 1 and the remaining 43,270 with 0. \n",
    "We draw half the marbles (21,869) from the urn to receive the treatment, and\n",
    "the half remaining in the urn receive the placebo. \n",
    "\n",
    "The results of the experiment are simply the count of the number of marbles marked 1 that were randomly drawn from the urn. \n",
    "\n",
    "We can simulate this process to get a sense of how likely it would be under these assumptios to draw only 117 marbles marked 1 from the urn. Since we draw half of the marbles from the urn, we would expect about half of the 468, or 234, to be drawn. The simulation study gives us a sense of the variation that might result from the random assignment process. That is, we can approximate the chance that there would be 117 or fewer cases of the virus in the treatment group. This calculation is based on the above assumptions, namely that the vaccine has no impact.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Simulation study\n",
    "\n",
    "We saw earlier that we can generate random outcomes according to the hypergeometric probability distribution, rather than having to program the urn sampling from scratch.   \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulations_fast = np.random.hypergeometric(ngood=468, nbad=43270, nsample=21869, size=500000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "234.000648"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(simulations_fast)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Count of Cases in the Treatment Group')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAEtCAYAAAC75j/vAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABXiUlEQVR4nO3deXwN1/vA8U8WQRNb0KgQgt4gIS4htSRIxB5LEZQQjbV2YldLqSVaFLWVUFvUToMitooWVUoVbUlssQchElnn94ffna8ri5tKxI3n/Xp5tTPz3Jlz5t7MMzPnzBkTRVEUhBBCCCNgmtMFEEIIIQwlSUsIIYTRkKQlhBDCaEjSEkIIYTQkaQkhhDAakrSEEEIYDUla74Bff/0VBwcHXF1dSUhIyOniZAlfX188PDze+HajoqKIjY3NsvWlpKQQHByMj48P1atXp1q1arRo0YI5c+bw5MmTLNtOVnt5/48ePRoHB4c3Wobjx4/j4OCg969y5crUqlWLLl26sH379lSf2bJlCw4ODhw/fjzT27t+/bpBcQ4ODowePTrd6azwclly6u8hJ5jndAFE9gsJCeG9997j0aNHHDhwgKZNm+Z0kV5b3759iYuLe6PbPHz4MAEBAWzdupX33nsvS9Y5cuRIdu3aRbNmzfD29sbU1JRz586xbNkyfvrpJ9avX0+RIkWyZFvZqWPHjtSuXTtHtu3l5YWXlxcASUlJREVFERoaysiRIzl16hSTJ09WY2vWrElgYCDly5fP1DYmTJhAREQEq1evfmVsYGAgdnZ2matEJmzevJnJkydz9uxZdV5O/D3kFElauVxCQgJ79+6ldevWhISEsHXr1lyRtOrWrfvGt3n27FkeP36cZes7deoUP/74I6NHj6ZHjx56y9zd3RkyZAjLli1jxIgRWbbN7KLVatFqtTmybQcHB1q3bq03r2fPnowaNYr169fj6upK8+bNAShdujSlS5fO9DbCwsKwtbU1KPblsmS13377jfj4eL15OfH3kFPk9mAud/jwYR4/foyrqyv16tXjyJEj3Lt3L6eLJYDTp08DaR9wmjVrho2NDX/88ccbLlXuYGpqysSJEylUqBDfffddThdHZCFJWrncjz/+iImJCTVr1sTLy4vk5ORU9/o9PDwYN24cGzduxNPTk2rVqtGpUyeOHTv2n+PGjx/P2LFjqVKlCu7u7jx48ACAkydP4ufnp56Zd+vWjd9++0397M8//4yDgwODBw/WW+fnn3+Og4MDP//8M5D6Hr6vry99+vQhNDSUVq1aUaVKFVq0aMHhw4eJiYlhwoQJ1KpVi9q1azNhwgSePXumflZRFIKDg2nfvj1arZYqVarQtGlTli5dim6Us9GjR7NgwQIAPD098fX1VT9/6dIl+vfvj4uLC87OznTq1IkjR4688ruxtLQEYMOGDaSkpKRaHhoaytq1a/Xm/frrr/Ts2RNXV1ccHR1xc3NjwoQJeleAo0ePpmXLlvz+++907NiRqlWr4unpydatW0lMTOTrr7+mbt261KpViyFDhvDw4UO9/ejn58eBAwdo3rw5VatWpU2bNuzZsyfDurzcpjV69GiaNm3K2bNn6dq1K87OztSpU4epU6fq7XuA8PBw+vXrh4uLC66urkydOpUNGzbg4ODAjRs3Xrkf02NlZUXDhg05f/489+/fB9Ju09qzZw/t2rVDq9VSo0YNevTowe+//64ud3BwIDIykhMnTuDg4MCWLVu4ceMGDg4OrFy5ks6dO+Pk5ISfn58an1Yb1uLFi3Fzc8PZ2Zlu3brp3d7L6HMvzvf19WXr1q1pzn+5Tevvv//ms88+w8XFhapVq+Lj40NoaKhejK+vL/7+/vz88898/PHHVKlShQYNGjB//vw0f5NvA0lauVhMTAyHDh2iWrVqFCtWjPr162NhYaH+6F/0yy+/8MUXX9CkSRMGDx7MgwcP6NmzJydOnPhPcTt37uTixYuMGzcOHx8frK2t2b9/P76+vty6dYt+/frRr18/bt26hZ+fH/v37wee3xZr27YtP/30k3rgP3r0KBs2bKBTp064u7unW9+//vqLsWPH0rhxYwICAnj48CFDhgyhd+/eREZGMnToUOrUqcMPP/zAsmXL1M/NnTuXSZMmUaFCBcaMGcOwYcPImzcvX3/9Ndu2bQOet9no2k3GjBlD3759gecHho4dO3Lp0iX69OnD0KFDSUpKonfv3uzatSvD76dx48YUKlSI1atX06hRI2bOnMnPP/+sdvSwsLDQiw8LC+PTTz8lLi6OQYMGMW7cOKpWrcoPP/zA9OnT9WLv3btH3759qVGjBqNGjcLc3JyxY8fSp08fjh07xmeffUbLli3ZvXs3gYGBep+9fPkygwYNombNmgQEBGBqasqgQYP48ccfM6zPyx48eIC/vz/lypVj3LhxVK9endWrVzNv3jw15ubNm3zyySecPn2aTz/9FH9/f/bt28fXX3+dqW2l58MPPwSef09pOXHiBEOHDqV48eKMGjWKAQMGcO3aNXr06KF2dggMDKRIkSKUK1eOwMBAatasqX7+m2++wcbGhrFjx+Lt7Z1uOfbs2cOKFSvo1KkT/fv3Jzw8nG7duvHvv/9mqj59+/bFxcVFLVfHjh3TjDt79iwdO3bk7Nmz9OjRg2HDhpGYmEj//v1TnQj9888/DBkyBFdXV8aPH0/p0qVZsGABwcHBmSrbG6OIXGvTpk2KRqNRli9frs7r3bu3otFolDNnzqjzGjZsqGg0GmXfvn3qvKioKMXFxUXx8fH5T3EVK1ZUrl69qs5LTExU3N3dlfr16ytPnjxR50dHRytubm6Km5ubkpCQoCiKojx69EipW7eu4uXlpTx48EBp0KCB0qhRI+Xp06fq57p27ao0bNhQb1qj0SgHDhxQ561Zs0bRaDR6ZUtJSVHc3d2Vjh07KoqiKAkJCUr16tWVoUOH6u27J0+eKE5OTkqfPn3UefPmzVM0Go1y/fp1ve2+XLbExETlk08+UerUqaPEx8crGfnjjz8UT09PRaPRqP8cHR2VPn366H1HiqIo/v7+SsOGDVOt08fHR9Fqter0qFGjFI1Go6xevVqdd+jQIUWj0aT6fKdOnZR69eql2o8rVqxQ58XFxSleXl5KvXr1lOTkZDXuxf2v2+bL06tWrdIra7NmzfS2N2bMGKVy5crKpUuX1Hm3b99WqlWrlmpfv+zYsWOKRqNR5s2bl27Mhg0bFI1Go4SEhCiKoiibN29WNBqNcuzYMUVRFGXixImKVqtVUlJS1M9cvHhRady4sbJ79251XsOGDZWuXbuq09evX1c0Go3i5eWl/m51NBqNMmrUKL3pSpUqKRcvXlTnXblyRXF0dFQGDBiQ7ufSm//yvlaU1N9Hhw4dlGrVqim3bt1S5z179kxp27atUrVqVSUqKkr9nEajUfbv368XV7NmTfVv5G0jV1q5mO7MWHeF8OL/v3y1Va5cORo1aqROW1tb07p1a86cOUNUVFSm4+zs7PR6UJ0/f57bt2/TpUsXrKys1PkFCxaka9eu3Llzh3PnzgFQqFAhJk+ezNWrV+nQoQO3b99m5syZr+yxlzdvXtzc3NRpe3t74PntPB0TExNsbW3Vdr08efKoV48vevjwIVZWVhl2b3/48CEnTpygfv36PHv2jAcPHvDgwQMeP36Ml5cX9+/f588//8ywzM7Ozvz0008sWbKEjh07UqpUKRITEzl48CAdO3bUu7pZsmQJmzdv1rsCy6icL37vZcuWBcDNzU3v86VKlUrVxlmgQAE++eQTdTpfvnx07tyZu3fvqt+RoZo1a6Y3XbFiRfV3oigK+/fvx83NTa83n42NDa1atcrUdtKTmJgIPP/e01KiRAmePn3K1KlTuXz5MvD8ttuePXsM6rD00UcfkSdPnlfGubm56d0+LVOmDG5uboSFhZGcnGxIVQx2//59zpw5Q+vWrSlRooQ6P2/evPj7+/Ps2TN++eUXdX7+/Plp0KCBXpy9vb16S/VtI70Hc6m7d+9y4sQJypYti4mJido2ULFiRUxMTNi5cydjxoxRD2AVKlRItY4yZcqgKAqRkZEULVo0U3G6/+rotq9LJC8qV64c8PxWka4HmqenJ40bN2bv3r107tyZ6tWrv7LOhQsXxtz8fz9pMzOzNMtiZmamtlXB88R16NAh9u/fT0REBFevXiU6OhpAL+5luttHq1evTrcr9K1bt15ZbnNzcxo0aKAeOMLDw1m3bh2rV69m6tSpeHl5kS9fPszMzLh+/TrffPMNly5d4tq1a9y5cyfd9b5Yb0P3BTw/4Xj51mSZMmUAiIyMpGrVqq+sk461tbXetIWFhXqQfvToEY8ePVIT6ot0v4nX9ejRI4B0Hxvo2rUrYWFhrFmzhjVr1lCqVCkaNmxI+/btqVix4ivX/3L90pNWfezs7Dhw4AAPHjygePHiBq3HEJGRkUDaf2u6k4ObN2+q8woXLoypqf71i4WFxVvbpiVJK5fatWsXycnJXLlyRe9KQyc6OprQ0FC1K3BaZ4u6g4vugJeZuBf/HzI++OuWvbju2NhYzp8/D8CRI0eIjY195ZXWiwnrRemdZeu2PWLECEJCQqhRowZarZaOHTtSs2ZNunfvnuH2dPXu0qWL3tXni9JK8joLFizAxsaGDh066M0vV64c48ePJzExkfXr13Pp0iWcnJxYv349EydOxN7eHhcXFxo3boyzszOrV69Os70prf2R0b7QSes71h3AXv5eX+Xlg+GLkpKSgNRtd/D8bD8rXLhwARMTk3QffLaysmLNmjX88ccfhIaG8vPPP7N69WrWrl1LYGBghu1UkPn98SJD9ul/uQrL6G9Nt80Xv+OMvqO3kSStXErXa3DGjBl6t+MALl68yPz589m6dauatK5du5ZqHVevXsXMzIxSpUqp8wyNe5nuGZfw8PBUyyIiIgD0bmXMnj2byMhIRo4cyaxZs5g9ezbjx4/PqMr/ycmTJwkJCeGzzz7T67GYlJTEo0ePMnymR1cnMzMz6tSpo7fs0qVL3Lhxg/z586f7eV0nj/bt26eZTDQaDfD89k18fDwzZszA1dWVoKAgvYT0zTffvLqimXDjxg0URdEr05UrV4D/XXFlhaJFi/Lee++p637R1atXX3v9MTExhIWFodVq070iioiI4MmTJ1SrVo1q1aoREBDApUuX6NKlCytWrHhl0jKU7urnRVevXqVAgQLqVaCpqWmqEWv+yy26zP6tGRvjSrHCIFeuXOHcuXPUqlWLNm3a0KhRI71/ffr0oXjx4hw9elS9vfTnn3/qPRN0//59duzYwUcffUShQoXU+YbGvczR0ZHixYsTHBxMTEyMOj8mJoZ169ZRvHhxnJycAPj9999Zu3YtPj4++Pv7065dO9asWcPJkyezaA/9j+720ctXRBs2bCAuLk69GoD/nZHqzmTff/99nJyc2Lp1q95tusTERMaOHcugQYP0Pv8yb29vrl+/zuLFi1Mti4+PZ9u2bZQtW5Zy5crx7Nkz4uLiKFu2rF7CunDhgtpzM6NtZcb9+/fZvXu3Oh0XF0dwcDBly5bN0qGaTE1N8fDw4Oeff9Yblig6OpqQkJDXWreiKEybNo3Y2Fh69uyZbtzUqVP57LPPePr0qTqvXLlyFCxYUO8KxNTU9LVulx05ckTvN/LPP/8QFhaGh4eHenJQrFgxLl68qHellFYPVF250iuP7m9px44d3L59W52fkJDAihUrsLCwMOqHkeVKKxfS3Spq3759msvz5MlDu3btWLx4sfrMloWFBb169aJ79+7ky5ePdevWkZKSwsiRI/U+a2hcWtv8/PPPGTJkCO3atVPLtmnTJu7evcu8efMwNTUlPj6ecePGYW1tTUBAAAABAQGEhoYybtw4tm/fTr58+V5r/7xIq9ViZWXF9OnTuXnzJgULFuT48ePs2rWLvHnz6h3MdGfry5Ytw93dHU9PT8aPH0/37t1p164dnTt3pnDhwuzcuZMzZ84wfPjwDIdg6tOnD8ePH2fu3LkcPnwYT09PrK2tuXXrFj/++CO3b98mKCgIExMTChUqhLOzM1u2bMHKygp7e3v+/fdfNm7cqB7Enj59muGJg6Hy5MnDmDFj+Ouvv3j//ffZvHkzd+7cSTO5vq7Bgwdz+PBhOnbsiK+vLxYWFqxfv1597syQ25l///23+jtOTk7m/v37hIaGcubMGbp165bm7XGdHj160KtXL7p06UKbNm3ImzcvoaGhXLt2jZkzZ6px1tbWXLx4kXXr1lGrVq1M/wYtLCz45JNP8PX1JS4ujpUrV1KwYEGGDBmixrRs2ZKgoCAGDBhAgwYN+Ouvv9i9e3eqq0Td9Lx583B1dU1z+Czd77J9+/Z07twZS0tLduzYwV9//cX48eMpWLBgpsr/NpGklQuFhIRQoEABGjdunG6Mj48PS5cuVXsR6gZqXbhwIU+ePMHFxYXhw4enaow2NC4tTZo0ISgoiIULF/Ltt99ibm6Os7MzX375pfrsyfz584mIiGDWrFnqH1aRIkUYMWIE48aN45tvvmHUqFH/ddekUqxYMZYuXcpXX33FwoULsbCwwN7entmzZ3P27FlWrVrF/fv3KVasGC1atGDv3r1s2bKFEydO4OnpiVarJTg4mPnz57NixQqSkpKwt7dnxowZtG3bNsNt58uXj1WrVhEcHMzu3btZtmwZT58+xdramjp16tCnTx+9xvRvvvmG6dOns3nzZhISErC1taV3796UL1+egQMHcuzYMZo0afLa++T9999n7NixzJw5k3v37uHo6MiKFSv0nk/KKnZ2dqxZs4aZM2eyZMkS8ubNS5s2bTAzM2P58uVptne9bN++fezbtw94nnDff/997O3tmTNnjnr7Oz316tVj0aJFLFmyhIULFxIfH8+HH37I7NmzadGihRo3cOBAJk6cyLRp0+jfv3+mbxt27NgRExMTFi9eTHx8PK6urowePZqSJUuqMYMHDyYpKYmdO3cSFhaGs7Mz33//vXryptO5c2eOHTvGsmXL+PPPP9NMWrrf5bx58wgKCiIlJYWKFSvy7bffptv+aixMlIxa7cQ7wcPDA1tb21cOBmponDBevr6+REZGcuDAgTeyvaioKKytrVNdUU2ZMoXg4GDOnDljUJdy8e6QNi0hRI4ZPHgwLVq00GufiYuL4+DBg1SsWFESlkhFbg8KIXJM69atGT9+PL1798bT05P4+Hi1A8GLrxQRQkeSlhAix3To0IG8efOyatUqZs2ahampKU5OTqxcuZJatWrldPHEW0jatIQQQhgNadMSQghhNOT2YDZSFIV36TpW1wHsXamz1Df3e9fq/DbU18Qk4+fzJGlls6iomFcH5RKFCj0fsig6Oi6HS/JmSH1zv3etzm9DfYsWtSKjZ8rl9qAQQgijIUlLCCGE0ZCkJYQQwmhI0hJCCGE0JGkJIYQwGpK0hBBCGA1JWkIIIYyGPKclxFvO0jIv5uaGn18mJaXw9Gl8NpZIiJwjSUuIt5y5uSkJyQoRkdGvjLW3LYRFJhKcEMZGkpYQRiAiMpqxi46+Mm5av7o42BXO/gIJkUPklEwIIYTRkKQlhBDCaEjSEkIIYTQkaQkhhDAakrSEEEIYDUlaQgghjIYkLSGEEEZDkpYQQgijIUlLCCGE0ZCkJYQQwmhI0hJCCGE0JGkJIYQwGpK0hBBCGA1JWkIIIYyGJC0hhBBGQ5KWEEIIoyFJSwghhNGQpCWEEMJoSNISQghhNCRpCSGEMBqStIQQQhgNSVpCCCGMRo4nrd9++43OnTvj7OxMvXr1mDJlCk+fPlWXh4WF0a5dO5ydnfHw8CAoKCjVOv788098fX3RarXUq1eP2bNnk5iYqBdz5coV+vbti4uLC66urkycOJGYmBi9mPv37zN8+HBcXV2pUaMGw4YN4969e9lTcSGMgKVlXgoVyk+hQvkxNzfD3NxMnX75n6Vl3pwurngHmOfkxv/44w969OiBh4cHixYt4urVq8yePZsHDx4wZ84cTp06Rd++fWnWrBmDBw/m999/JzAwEEVR8Pf3B+Dq1av4+fmh1WqZO3culy9fZs6cOcTExDBhwgQAoqOj6d69O8WLF2fmzJlERUUxa9Ysbt++zZIlSwBISkrC39+f2NhYJk2aRFJSEl9//TU9e/Zk8+bNmJvn6K4SIkeYm5uSkKwQERmdYZy9bSEszHP8HFi8A3L0SPzVV19RrVo1vvnmG0xMTKhTpw4pKSmsWLGCuLg45s2bR+XKlZk1axYA7u7uJCUlsXjxYnx9fbGwsGDp0qUUKFCAhQsXYmFhQf369cmXLx9Tp06lT58+2NjYsHbtWh4/fsy2bdsoUqQIADY2NvTu3ZszZ87g7OzMzp07uXjxIrt27aJ8+fIAVKpUiZYtW7J3716aN2+eY/tJiJwUERnN2EVHM4yZ1q8uDnaF30yBxDstx06NHjx4wMmTJ+ncuTMmJibq/C5duhAaGoqpqSknT56kcePGep9r0qQJjx8/5tSpUwAcPXqUhg0bYmFhocY0bdqU5ORkwsLC1JiaNWuqCQugXr16WFpacvjwYTWmQoUKasIC1GldjBBCiJyVY1da//zzD4qiUKhQIYYMGcKhQ4cwMzOjZcuWjBkzhhs3bpCYmIi9vb3e58qUKQNAREQEzs7O3Lp1K1WMtbU1VlZWREREABAeHk6rVq30YszMzChVqpRezMvrAbCzs1Nj/otChfL/588aG3NzM+DdqfObqq9uO5mJz6oyZWbbWbndt4X8pt+8F65h0pSjV1oAo0ePpkiRIixatIiBAweyfft2Jk2axJMnTwCwsrLS+5ylpSUAMTEx6cbo4nQdLZ48eZIlMUIIIXJWjl1p6Xr3Va9enYkTJwJQu3ZtFEVh5syZ+Pj4AOjdOnyRqakpiqKkG6MoCqam/8vJWRWTWdHRcf/5s8ZGd3b2rtT5TdU3s2e9SUnJWVamzGw7K7f7tpDf9JtXtKhVhldbOXalpbticnd315tfr149FEXhzz//BEh1laObLlCggHpllNaVUGxsLAUKFACeX4mlFfP06VN1HYbECCGEyFk5lrTKli0LQEJCgt583RVYqVKlMDMz49q1a3rLddP29vZYWlpiY2PD1atX9WKioqKIiYlR26js7e1TxSQnJ3Pjxo0MY3TbS6utSwghxJuXY0mrfPny2NrasmvXLr35Bw8exNzcHK1Wi4uLC3v37lVvAwLs2bOHAgUK4OTkBEDdunU5ePCgXvLbs2cPZmZm1KpVS405fvw4jx49UmPCwsKIjY2lTp06wPMrvH///Zfw8HA15tKlS4SHh6sxQgghclaOJS0TExMCAgI4efIkAQEB/PLLLyxdupRFixbh6+uLtbU1/fr149SpUwwdOpTDhw8zd+5cli9fTp8+fcif//m91549e3Lv3j169+7NwYMHWbFiBdOnT8fHx4eSJUsC8Mknn2BhYYGfnx/79u1j48aNjBgxAnd3d6pXrw5A8+bNKVOmDD179mTnzp2EhITQq1cvPvzwQ5o1a5ZTu0kIIcQLTJQXL2NyQGhoKN9++y2XLl2iaNGidOzYkT59+qidH/bt28e8efOIiIjAxsaGLl268Omnn+qt4+TJkwQGBnLhwgWKFClCmzZtGDhwIHny5FFj/vnnH6ZNm8bp06extLSkUaNGjBw5Uq+96tatW3z55ZccPXoUCwsL6taty+jRo3n//ff/U90UReH+/Xen5+Hb0Ij7Jr3Jjhh/X3v0ygd84X8P+WZlRwxDtp3V231byG/6zSta1ApT0/R7YmSq92BMTAz//vsvWq0WeJ4sVq1ahZmZGV26dMHFxSXTBWzUqBGNGjVKd7mXlxdeXl4ZrsPFxYUNGzZkGKPRaFi5cmWGMR988AELFizIMEYIIUTOMThpXbp0iW7dulG0aFF+/PFHrl+/To8ePVAUhTx58rBv3z6+++47ateunZ3lFUII8Q4zuE1r7ty5AIwYMQKAjRs3kpSUxOrVq/nll1+oVKkSixYtypZCCiGEEJCJpPXbb7/h5+enPld14MABypQpg1arJX/+/LRp04Zz585lW0GFEEIIg5NWfHy8OuBsZGQkly5dws3NTS/GzCxzY6QJIYQQmWFwm5adnR2nTp2iQ4cObN26FRMTEzw9PYHnveR++ukndTBbIUTO+KCYpcED1yYlpfD0afwbKJUQWcfgpNW5c2cmT57MuXPnCA8P58MPP+Sjjz7in3/+YdSoUVy8eJEZM2ZkZ1mFEK+QL685sfFJ8tJGkWtlKmlZWloSEhKCVqulf//+6rJnz54xZcoUWrdunS2FFEIYTl7aKHKzTD2n1apVq1TvpdJoNOzevTtLCyWEEEKkJdOvJomPj+fRo0ckJyenuVw3dJIQQgiR1QxOWo8ePWLy5Mns27cv3YQFcOHChSwpmBBCCPEyg5PWjBkz2L17N25ublSqVAkLC4vsLJcQQgiRisFJ68CBA3To0IEpU6ZkZ3mEEG+AoV3jzc3l2UvxdjE4aSUlJVGlSpXsLIsQ4g0xtGt8ZXvrN1QiIQxjcNKqWbMmx48fx8fHJzvLI4R4QwzpGh88tfkbKo0QhjE4aY0dO5Zu3boRGBhI06ZNsba2Vt959SLpPSiEECK7GJy0vL29SUlJISgoiBUrVqQbJ70HhRBCZBeDk1avXr0wMUn/bZJCCCFEdjM4aQ0cODA7yyGEEEK8UqZHxDhy5AihoaHcvHmTPHnyULJkSRo0aEC9evWyo3xCCCGEyuCklZKSQkBAALt370ZRFAoWLEhKSgoxMTGsXbuWxo0bM3fuXLmFKIQQItsY/G6CZcuWsWvXLjp37kxYWBgnTpzg5MmThIWF0bVrV/bs2cP333+fnWUVQgjxjjM4aW3ZsoVGjRoxYcIEihUrps4vVqwY48aNw8vLi02bNmVLIYUQQgjIRNKKjIykbt266S6vXbs2169fz5JCCSGEEGkxOGkVKVKEK1eupLv8ypUrFChQICvKJIQQQqTJ4KTl4eFBcHAwBw4cSLVs//79rF+/Hg8PjywtnBBCCPEig3sPDhkyhF9//ZX+/ftTvnx57O3tAQgPDyc8PBxbW1uGDBmSXeUUQgghDL/SKly4MBs3bsTf3x9FUfj55585fPgwKSkp9OjRg82bN2NtLSNCCyGEyD6Zeri4YMGCBAQEEBAQkF3lEUIIIdKVbtK6efMm1tbW5MuXT502hIzyLoQQIrukm7Q8PT0JDAzE29sbeN4Rw5DRLmSUdyGEENkl3aTVv39/HBwc9KZliCYhhBA5Kd2kNWDAAL1pQ0Z5T0hIeP0SCSGEEOkwuPegp6cn+/fvT3d5SEgIbm5uWVIoIYQQIi3pXmk9ePCAy5cvq9ORkZH8+eefFCxYMFVsSkoK+/btkystIYQQ2SrdpJU3b16GDx/OvXv3ADAxMWHJkiUsWbIkzXhFUWjevHn2lFIIIYQgg6RlaWnJokWL+Oeff1AUhbFjx+Lj44NWq00Va2pqirW1NbVr187WwgohhHi3ZfhwsaOjI46OjsDz57QaN26MRqN5IwUTIreztMyLufmrm5XNzc3eQGmEMA4Gj4jxcm/CtJw/f57KlSu/VoGEeFeYm5uSkKwQERmdYVxlexkeTQgdg5NWYmIiS5cuZe/evcTGxpKSkqIuS05O5unTp8TExMjDxUJkQkRkNGMXHc0wJniqtBULoWNwl/e5c+cyf/58oqOjyZ8/P5GRkXzwwQeYm5tz+/ZtEhMTGTduXHaWVQghxDvO4KT1008/UatWLQ4cOMB3330HwIQJE9izZw9LliwhKSmJPHnyZFtBhRBCCIOT1p07d2jcuDGmpqbY2NhQtGhRTp8+DUD9+vVp27YtGzZsyLaCCiGEEAYnrXz58uldSdnZ2fHPP/+o01WrVuX69etZWzohhBDiBQYnrUqVKvHzzz+r0+XKlVOvtOD5lZgMqCuEECI7GZy0unTpwv79+/nkk0+IiYmhRYsWnD9/njFjxvDdd9+xcuVKqlSpkp1lFUII8Y4zuMt706ZNmTJlCitWrCB//vzUqVOHXr16qZ0ySpYsyZgxY7KtoEIIIYTBV1oAHTp0YNeuXZiZPX9Cf/jw4Rw4cICtW7fy008/Ub58+f9ckAEDBuDl5aU3LywsjHbt2uHs7IyHhwdBQUGpPvfnn3/i6+uLVqulXr16zJ49m8TERL2YK1eu0LdvX1xcXHB1dWXixInExMToxdy/f5/hw4fj6upKjRo1GDZsmDruohBCiLeDwVda6SlZsiQlS5Z8rXVs376dffv2YWdnp847deoUffv2pVmzZgwePJjff/+dwMBAFEXB398fgKtXr+Ln54dWq2Xu3LlcvnyZOXPmEBMTw4QJEwCIjo6me/fuFC9enJkzZxIVFcWsWbO4ffu2OvhvUlIS/v7+xMbGMmnSJJKSkvj666/p2bMnmzdvxtz8tXeTEEKILJDu0djT0zPTKzMxMSE0NDRTn7lz5w5ffvklJUqU0Js/b948KleuzKxZswBwd3cnKSmJxYsX4+vri4WFBUuXLqVAgQIsXLgQCwsL6tevT758+Zg6dSp9+vTBxsaGtWvX8vjxY7Zt20aRIkUAsLGxoXfv3pw5cwZnZ2d27tzJxYsX2bVrl3q1WKlSJVq2bMnevXtl9HohhHhLpHt7UHcFlZl/H3zwQaYLMH78eOrWras3Qnx8fDwnT56kcePGerFNmjTh8ePHnDp1CoCjR4/SsGFDLCws1JimTZuSnJxMWFiYGlOzZk01YQHUq1cPS0tLDh8+rMZUqFBB7/ambloXI4QQIuele6W1evXqbN/4xo0b+euvvwgJCSEwMFCdf/36dRITE7G3t9eLL1OmDAARERE4Oztz69atVDHW1tZYWVkREREBQHh4OK1atdKLMTMzo1SpUnoxL68Hnj+LposRQgiR83KssSYyMpLp06czffp0rK31R7F+8uQJAFZWVnrzLS0tAYiJiUk3Rhen62jx5MkTg2IqVKiQZszVq1czWzU9hQrlf63PGxPdKzTelTq/bn1z2ytHzM3Nct13L7/pN+9Vj/sanLQMbePav3//K2N0L5WsX78+TZo0SXM5kO7DyqamphnGKIqCqen/7nxmVYwQQoicZXDSSquHYEpKCvfv3+fq1auULVuWunXrGrSutWvX8vfff/Pjjz+SlJQE/C9RJSUlUaBAAYBU3dJ10wUKFFCvnl6OAYiNjVXXYWVllWbM06dPsbW1fWVMWldpmREdHfdanzcmurOzd6XOr1vf3Hb2npSUnOu+e/lNv3lFi1pleLVlcNLKqI3r3Llz9OzZk1q1ahm0rj179vDw4UPq1auXapmjoyOTJk3CzMyMa9eu6S3TTdvb22NpaYmNjU2q23dRUVHExMSobVT29vapYpKTk7lx44Z6lWdvb683juKL23N2djaoTkIIIbJfltz7cnJyomvXrnz77bcGxU+ePJlNmzbp/WvYsCElSpRg06ZNNG3aFBcXF/bu3ategcHzZFegQAGcnJwAqFu3LgcPHiQhIUEvxszMTE2gdevW5fjx4zx69EiNCQsLIzY2ljp16gDPexP++++/hIeHqzGXLl0iPDxcjRFCCJHzsqwjRrFixbhy5YpBseXKlUs1r3DhwlhYWKjjF/br148ePXowdOhQ2rZty+nTp1m+fDnDhw8nf/7nl7A9e/Zk586d9O7dm+7du3PlyhVmz56Nj4+Pejvzk08+Yc2aNfj5+dG/f38ePXrErFmzcHd3p3r16gA0b96cxYsX07NnT4YPH46iKHz99dd8+OGHNGvWLAv2jhBCiKyQJVda9+7dIzg4+LVHxnhR7dq1mT9/PpcvX6Z///78+OOPjBw5kl69eqkx5cuXJygoiNjYWAYNGsSKFSvo0aOH3huUra2tWbVqFYULFyYgIIA5c+bQtGlT5syZo8ZYWFiwYsUKKleuzPjx45kyZQparZbly5fLaBhCCPEWee3egwkJCTx48IDk5GQmTpz4nwsyY8aMVPO8vLxSjUf4MhcXl1e+fFKj0bBy5coMYz744AMWLFjwynIKIYTIOa/VexCeP6jr6upKy5YtadCgQVaVSwghhEglS3oPCiGEEG+CPDkrhBDCaBh8pfXo0SMCAwM5evQo9+7d0+uKrmNiYsL58+eztIBCiLffB8UsDR7GKSkphadP499AqURuZHDSmjx5Mrt376Z69eq4urqqL4IUQuiztMyLufmrb2LkprEH8+U1JzY+iYjI6Azj7G0LYWHAvhEiPQYnrV9++YWuXbsyfvz47CyPEEbP3NyUhGTllQfwyvbWGS43NhGR0YxddDTDmGn96uJgV/jNFEjkSgYnrTx58qT5ULAQIjVDDuDBU+XlokJklsHX6W3btmX79u3qALdCCCHEm2bwldbgwYPp06cPTZo0wd3dnaJFi6aKMTExoX///llaQCGEEELH4KQVEhLCr7/+SkpKCsHBwWnGSNISQgiRnQxOWgsWLMDOzo4xY8Zgb28vvQfFO+flXoHpveU1N/UKFOJtY3DSunfvHqNHj6Z+/frZWR4h3lrvaq9AId4mBietSpUqERkZmZ1lEeKtJ70ChchZBvceHDlyJBs3bmTdunXcvXuXlJSU7CyXEEIIkYrBV1q6145MmTKFKVOmpBkjwzgJIYTITgYnLUdHR/U190IIIUROMDhppfWSRiGEEOJNkpErhRBCGI10r7QqVapEYGAg3t7eAFSsWBETE5MMVyZtWkIIIbJTukmrTZs22NnZ6U2/KmkJIYQQ2SndpDV9+nS9aWnTEkIIkdNeq00rOjqa+Hh5A6kQQog3I8OklZiYyPr16xkzZoze/JMnT9KiRQs++ugjtFotPXv25Nq1a9laUCGEECLdpJWQkED37t2ZNGkSISEh6nu0rly5gr+/P+Hh4bi5ueHn50dERASdOnXi/v37b6zgQggh3j3pJq3vv/+e06dPM2LECH777TfMzZ83f82fP5/4+HhatGjB0qVLGTlyJJs3b8bMzIzFixe/sYILIYR496SbtHbv3k2TJk3w9/cnX758wPOrrwMHDmBiYoK/v78aW7hwYT7++GMOHTqU7QUWQgjx7ko3aV29ehUXFxe9eX/88QdxcXEUL16cSpUq6S2zs7Pj7t272VNKIYQQggySVkpKSqoXPf76668A1KlTJ1X8kydPyJ8/f6r5QgghRFZJN2nZ2dlx4cIFvXmhoaGYmJjQoEGDVPFhYWF6DyMLIYQQWS3dpNWiRQu2b99OaGgocXFxrFy5kn///ZeiRYvi4eGhF7tjxw6OHj2Kp6dnthdYCCHEuyvdETH8/Pw4cuQIAwYMwMTEBEVRyJMnD19++SUWFhYA7Nu3jzVr1nDixAns7e3x8/N7U+UWQgjxDko3aVlYWLBy5Up27drFH3/8gaWlJa1ataJChQpqzLlz5zh16hStWrVi9OjRai9DIYRIywfFLDE3N6NQoVe3fyclpfD0qYy4I/Rl+D4tMzMzvL291ZHeX9a3b18GDx6Mqam84UQI8Wr58poTG59ERGR0hnH2toWwMJfjikjN4JdApkV6CwohMisiMpqxi45mGDOtX10c7Aq/mQIJoyKnMkIIIYyGJC0hhBBGQ5KWEEIIo5Fu0goODubKlStvsChCCCFExtJNWoGBgZw8eVKd9vT0ZP/+/W+kUEIIIURaMnxOKzQ0lGrVqpE/f34iIyO5efMmN2/ezHCFJUuWzPJCCiGEEJBB0mrfvj3Lly/n8OHDAJiYmDBt2jSmTZuW4QpfHq9QCCGEyCrpJq0RI0ZQs2ZN/v77bxISEvj222/x8vLCwcHhTZZPCCGEUGX4cHGDBg3UEd23bt1KmzZtZFBcIUS2k+GeRHoMHhHjwIEDACQnJ3Pu3DkiIyOxsLCgRIkSODk5ZVsBhRDvHhnuSaQnU8M4HTx4kMmTJ3Pnzh0URQGet3W9//77TJw4MdUrS4QQ4r+S4Z5EWgxOWidPnmTgwIEULVqUoUOHUr58eRRFITw8nHXr1jFo0CBWrVpF9erVs7O8Qggh3mEGX1fPnz8fW1tbQkJC6N27N56enjRq1IjevXsTEhKCra0tixYtytTGU1JSCA4OxtvbG61WS6NGjZg+fToxMTFqTFhYGO3atcPZ2RkPDw+CgoJSrefPP//E19cXrVZLvXr1mD17NomJiXoxV65coW/fvri4uODq6srEiRP1tgNw//59hg8fjqurKzVq1GDYsGHcu3cvU3USQgiRfQxOWmfPnqVDhw4UKFAg1TIrKyvat2/PmTNnMrXxZcuWMWXKFBo0aMC3335Ljx492LZtG4MHDwbg1KlT9O3bl3LlyjF//ny8vb0JDAxk+fLl6jquXr2Kn58fefPmZe7cuXz66aesWLGC6dOnqzHR0dF0796d+/fvM3PmTIYPH86uXbsYPny4GpOUlIS/vz9nz55l0qRJTJo0iVOnTtGzZ0+SkpIyVS8hhBDZ47VeTfIiExOTVFc3GVEUhWXLltGxY0c1edSpU4ciRYowdOhQLly4wLx586hcuTKzZs0CwN3dnaSkJBYvXoyvry8WFhYsXbqUAgUKsHDhQiwsLKhfvz758uVj6tSp9OnTBxsbG9auXcvjx4/Ztm0bRYoUAcDGxobevXtz5swZnJ2d2blzJxcvXmTXrl2UL18egEqVKtGyZUv27t1L8+bNs2pXCSGE+I8MvtJydnZm06ZNxMbGploWExPDxo0bqVKlisEbfvr0Ka1ataJly5Z688uVKwfAv//+y8mTJ2ncuLHe8iZNmvD48WNOnToFwNGjR2nYsCEWFhZqTNOmTUlOTiYsLEyNqVmzppqwAOrVq4elpaX68PTRo0epUKGCmrAAdVoXI4QQImcZfKU1YMAAunXrRsuWLenatStly5YFUDti3Llzh8mTJxu8YSsrK8aPH59qfmhoKACVK1cmMTERe3t7veVlypQBICIiAmdnZ27dupUqxtraGisrKyIiItQytmrVSi/GzMyMUqVK6cW8vB4AOzs7NUYIIUTOMjhpubi4MH/+fL744gsCAwMxMTEBnt/mK168OHPmzOGjjz56rcKcOXOGpUuX0qhRI548eQI8T24vsrS0BJ5f3aUXo4vTdbR48uSJQTEVKlRIM+bq1av/uU6GPByZW5ibmwG5t866+om3i6EPIf/XdUPu/U2/7G2o7/+nlnRlqk3L09OTBg0a8Ndff3Hjxg0AbG1tcXR0xNz89ZrHfv/9d/r27UupUqWYOnWqenVjkk4NTE1N9Z4Ve5miKJia/u/uZ1bFCCGEyDmZzjRmZmZUrVqVqlWrZlkhdu3axejRoylbtizLli2jSJEi3L9/HyBVt3TddIECBdSrp5djAGJjY9WejlZWVmnGPH36FFtb21fGpHWVZqjo6Lj//Fljozs7y611flfOto1NUlJytv3mcvtv+mVvQ32LFrXK8Gorxy8hVqxYwbBhw6hWrRpr167l/fffB563JZmZmXHt2jW9eN20vb09lpaW2NjYpLp9FxUVRUxMjNpGZW9vnyomOTmZGzduZBij215abV1CCCHevBxNWhs3bmTGjBk0a9aMZcuW6T0DljdvXlxcXNi7d696GxBgz549FChQQB3vsG7duhw8eJCEhAS9GDMzM2rVqqXGHD9+nEePHqkxYWFhxMbGUqdOHeB5b8J///2X8PBwNebSpUuEh4erMUIIIXJWjiWtqKgovvzyS2xtbenSpQvnz5/njz/+UP89ePCAfv36cerUKYYOHcrhw4eZO3cuy5cvp0+fPuTP//wytmfPnty7d4/evXtz8OBB9cFiHx8f9YWUn3zyCRYWFvj5+bFv3z42btzIiBEjcHd3V4edat68OWXKlKFnz57s3LmTkJAQevXqxYcffkizZs1yajcJIYR4QZY9XJxZR44cIS4ujsjISLp06ZJqeWBgIK1bt2b+/PnMmzeP/v37Y2Njw8iRI/n000/VuPLlyxMUFERgYCCDBg2iSJEi9OjRg4EDB6ox1tbWrFq1imnTphEQEIClpSVNmzZl5MiRaoyFhQUrVqzgyy+/ZPz48VhYWFC3bl1Gjx792p1MhBBCZA2Dj8bdunWjX79+1K5dG3je+eGzzz5j9OjRVK5cOdMbbtOmDW3atHllnJeXF15eXhnGuLi4sGHDhgxjNBoNK1euzDDmgw8+YMGCBa8skxBCiJyRbtJyc3PD0dERR0dHKleuzIkTJ/Dx8VGXJyYmcuLECaKjM37fjRBCCJFV0k1a/v7+XLhwgb1797JkyRJMTEz44osv2LBhA5UqVaJ06dKYmJik+xyVEEIIkdXSTVp+fn7q/yckJFC1alUaNGiApaUlZ8+eZdOmTSiKQt++falUqRJOTk5UqVIl1XBJQgghRFYxqE1LNxitm5sb3t7eADx48IA6derQtWtXkpOT+euvv9i+fbskLSGEENkm3aTl4+NDpUqVcHR0pGLFioD+MEe6/69bt67aOUMIY2RpmRdz81c//SFjDwqR89JNWjVr1uTixYvs27ePBw8eYGJiwty5czl8+DAVK1akZMmS0qYlcgVzc1MSkhUiIjPuVFTZ3voNlUgIkZ50k9aIESPU/799+zYNGjTgww8/5NmzZ6xfv14dMHfUqFE4Ozvj5OSEk5OTjB4hjFJEZDRjFx3NMCZ4qrwIVIicZlCbVokSJYDno0bo2rRu3ryJh4cH7u7uxMXFsXnzZubOncv58+ezr7RCCCHeaQY/XFyyZEnee+89ddrKyoqSJUvy8ccfo9VqgbRHWxdCCCGyisFJ68CBA3rTBQsWTDXvdV7hIYQQQrxKjr+aRAghhDCUJC0hhBBGQ5KWEEIIoyFJSwghhNGQpCWEEMJoyNsNRa4lwzMJkftI0hK5lgzPJETuI0lL5GoyPJMQuYu0aQkhhDAakrSEEEIYDUlaQgghjIYkLSGEEEZDOmIIIYzWB8UsMTc3o1Ch/AbFJyWl8PRpfDaXSmQnSVpCCKOVL685sfFJr3ysAcDethAWBjy3J95ukrSEEEbNkMcaAKb1q4uDXeHsL5DIVnLaIYQQwmhI0hJCCGE0JGkJIYQwGpK0hBBCGA1JWkIIIYyGJC0hhBBGQ7q8CyHeCZl5EFkeQn57SdISQrwTDH0QWR5CfrtJ0hJCvDMMeRBZHkJ+u8nphBBCCKMhSUsIIYTRkKQlhBDCaEjSEkIIYTSkI4YwOpaWeTE3oHeXubnZGyiNEOJNkqQljI65uSkJycoruy5Xtrd+QyUSQrwpkrSEUTKk63Lw1OZvqDRCiDdF2rSEEEIYDUlaQgghjIYkLSGEEEZD2rTEW8HQHoEgvQJF9npxYF3dby29QXZlYN03T5KWeCsY2iMQpFegyF4ysO7bTZKWeGsY0iMQpFegyH4ysO7bS5LWS0JCQli0aBHXr1/H1taWPn360KZNm5wultGSB4GFEFlJktYLdu/eTUBAAN26dcPNzY3Q0FBGjRpFvnz5aNq0aU4XzyjJg8Ait5KXSuYMSVovmD17Ns2aNWPs2LEAuLm5ER0dzTfffCNJ6zXIg8AiNzK07cuhbBHy5ctj0B0HSW6vJknr/12/fp1r164xbNgwvflNmjRh9+7dXL9+ndKlS+dQ6d4+ad32S6unldz2E7mZoSdkCaRIx44sYqIoipLThXgbHD58mN69e7N9+3YqVqyozj9//jxt27blu+++w93dPVPrzO271sTExKA4RVFIScl4X5iaPl/Xq+IyE/uuxRlDGWXfvDpOZHxskSut//fkyRMArKys9OZbWloCEBMTk+l1GnpQz+1MTEwwMzNsXxgal5nYdy0uJ7f9tsfl5LYzU0aRPrkW/X+6q6KXE41uvqmp7CohhMhpciT+fwUKFABSX1E9ffpUb7kQQoicI0nr/9nb2wNw7do1vflXr17VWy6EECLnSNL6f2XKlKFUqVL89NNPevP37t1L2bJlKVmyZA6VTAghhI50xHhB//79GTNmDIUKFaJBgwYcOHCA3bt3M2fOnJwumhBCCKTLeyrr168nKCiIW7duUbp0aXr37i3DOAkhxFtCkpYQQgijIW1aQgghjIYkLSGEEEZDkpYQQgijIUlLCCGE0ZCkJYQQwmhI0noNFy5cwNHRkdu3b+vN37NnD23atKFatWo0adKElStXkpKSohfz559/4uvri1arpV69esyePZvExMQ3WXyDpKSkEBwcjLe3N1qtlkaNGjF9+nS94a7CwsJo164dzs7OeHh4EBQUlGo9uam+u3fvpl27dmi1WurXr8+YMWOIiorSW09uqu+Lpk2bRuXKlVPNz031/fvvv/H390er1VK7dm1GjBjB/fv39dZjLPUFw+psVMcsRfwnly9fVtzc3BSNRqPcunVLnX/kyBFFo9EoI0eOVI4ePaosXLhQqVixovLdd9+pMVeuXFGqV6+u+Pv7K4cOHVKWL1+uODk5KZMnT86JqmRoyZIlSqVKlZSvvvpKOXr0qLJmzRqlVq1ayqeffqooiqL8/vvviqOjoxIQEKAcPnxYmT17tuLg4KAsW7ZMXUduqu/OnTsVjUajfP7558qRI0eULVu2KA0aNFBatmypxMfHK4qSu+r7ohMnTigVK1ZUKlWqpDc/N9X32rVrSo0aNZSuXbsqhw4dUrZt26a4ubkpPj4+6jqMqb6K8uo6G9sxS5JWJiUmJipr1qxRtFqtUqtWrVRJa9iwYUqjRo2U5ORkdd6IESOUhg0bqtNjx45V6tevrx7kFEVR1q5dq1SqVEm5ffv2m6mIAVJSUpSaNWsqkyZN0puvO3CfP39e6d69u9KhQwe95YGBgYqLi4tav9xU31atWim9evXSW/7HH38oGo1G2bdvn6Iouau+Ok+fPlU8PT0Vd3f3VEkrN9V35MiRipeXl/Ls2TN1+f79+xV3d3fl2rVriqIYT30VxbA6G9sxS24PZtLvv//OV199xaeffkpAQECq5fHx8eTPn1/vVSaFCxfm0aNH6vTRo0dp2LAhFhYW6rymTZuSnJxMWFhYtpY/M54+fUqrVq1o2bKl3vxy5coB8O+//3Ly5EkaN26st7xJkyY8fvyYU6dOAbmnvteuXaNOnTr4+PikuxxyV311Zs6cSbFixfj4449TrSc31Tc0NJT27duTN29edbmHhweHDx9W31xuLPUFw+psbMcsGXswk8qXL09oaChFixZly5YtqZZ36dKFXr16sXr1atq0acO5c+fYsmULrVu3BiAuLo5bt26lGjXe2toaKysrIiIi3kg9DGFlZcX48eNTzQ8NDQWgcuXKJCYmpqpLmTJlAIiIiMDZ2TnX1LdChQo0adIkw+W56futUKEC8PyAtX37drZu3UpISIhebG6rb0xMDCVKlGDChAns2rWLxMREPD09+fzzzylSpIhR1RcMq7OxHbPkSiuTihUrRtGiRdNd/tFHH/Hpp58ydepUXFxc8PPzo2rVqowZMwZI/w3J8Pwtyf/lDclv0pkzZ1i6dCmNGjUy6G3Puam+5cuXT7X82rVrzJw5E0dHR+rVq5fr6vvkyRPGjRvHoEGD0nw9T26qr+7deYGBgURHR/PNN98wfvx4jh49yqBBgwDjry+k/o6N7ZglV1pZbOLEiWzZsoUBAwbg6urKpUuXmDdvHoMHD2bhwoXpviEZnr8l+W1+Q/Lvv/9O3759KVWqFFOnTlXPsNKqCzx/23Nuqu/LLl++jL+/P+bm5sydOzdX1nfatGmUKFECPz+/ND+Tm+p7+fJlAGxsbJg7d65ap0KFCjFw4ECOHTumJm5jrC+k/R0b2zHr7d7DRubOnTts2LCBPn36MHDgQGrVqsUnn3zCV199xYEDBzh06JB6tpLW2UlsbOxb+4bkXbt20aNHDz744ANWrlxJkSJF0n3bs266QIECuaq+Lzp+/DidO3cG4Pvvv8fOzg4gV9X34MGD7Ny5kylTppCSkkJSUpLaDVr3/7mpvrq6uLu76x2g69atCzzvCm+s9YW062yMxyxJWlno5s2bKIpC9erV9ebXrFkTeN5xwdLSEhsbG/WNyDpRUVHExMS8lW9IXrFiBcOGDaNatWqsXbuW999/HwA7OzvMzMxSve1ZN21vb5+r6quza9cu/P39sbGx4YcfftC7bZib6rtnzx7i4+Np2bIljo6OODo6snDhQpKTk3F0dOTbb7/NVfUtXbo0JiYmJCQk6MUnJycDz680jLG+kH6djfGYJUkrC5UpUwYzMzN+++03vfmnT58GoFSpUsDzM7eDBw/q/XHs2bMHMzMzatWq9eYKbICNGzcyY8YMmjVrxrJly/TOqvLmzYuLiwt79+5VbyHA87oUKFAAJycnIPfUF+DIkSOMGDECrVZLcHAwNjY2qdaRW+o7YMAANm3apPfPx8cHMzMz9f8h99TX0tKSGjVqsG/fPr2HZg8cOACAi4sLYFz1hYzrbIzHLGnTykLW1tZ07dqVZcuWYWpqiqurKxEREcyfP5+KFSvSqFEjAHr27MnOnTvp3bs33bt358qVK8yePRsfHx9KliyZw7X4n6ioKL788ktsbW3p0qUL58+f11tuZ2dHv3796NGjB0OHDqVt27acPn2a5cuXM3z4cPLnzw/knvqWKlWKcePG8d5779G3b18uXbqkt/yDDz7AxsYm19TXzs5OPWjpHDp0CIAqVaqo83JTfYcOHYqfnx99+/bFz8+PGzdu8NVXX+Hl5aWOBGIs9QXD6mxsxyx5CeRr2LJlC2PGjOHw4cOUKFECeD5kysqVK/nhhx+IjIykRIkSeHh4MGjQIL3eNydPniQwMJALFy5QpEgR2rRpw8CBA8mTJ09OVSeVbdu2MWrUqHSXBwYG0rp1a/bt28e8efOIiIjAxsaGLl268Omnn+rF5ob6Dhw4kPnz56e7fPDgwXz22WdA7qiv7vt90fz581m0aFGqg19uqu/JkyeZPXs2586do0CBArRs2ZJhw4bpPbtlDPUFw+rs7e1tVMcsSVpCCCGMhrRpCSGEMBqStIQQQhgNSVpCCCGMhiQtIYQQRkOSlhBCCKMhSUsIIYTRkKT1lomJiSEoKIiPP/6YGjVqUK1aNdq3b88PP/yQ6vXXOSkqKorY2NgsW9/KlSupV68eVatW5auvvsowNiEhgQ0bNtCpUydcXV2pUqUK3t7eLF26NNUQPMZky5YtODg4cPz48Sxb5/Xr19X/P378OA4ODmm+Uicr6Mr/qn8eHh7Zsv2MvLgf3kaZKd+TJ09YuXIlPj4+fPTRR1StWpXGjRvzxRdfcPPmzWws5dtBntN6i4SHh9OvXz8iIyPx9vbGycmJ+Ph49u/fz8mTJ/H29mbWrFnpjqr+phw+fJiAgAC2bt2aasSE/+Lvv/+mVatWaoKuWLGi3ogLL7p37x79+/fn7NmzNG7cmOrVq2NmZsbRo0c5ePAgrq6uLFu2TO9ldcbi+vXrnDp1irp161KsWLHXXp+/vz/FixdnxowZwPOk1a1bN6ZPn57myxxfl678Og8fPmT69Om4uLjovTjT0tJSHWnhTZgwYQIRERGsXr36jW0zMzJTvkuXLvHZZ58RGRlJ06ZN0Wq15M2bl4sXL7J161YAli9fjlarze5i55w3+p5kka5nz54pTZo0UWrVqqVcuHAh1fJJkyYpGo1G+f7773OgdPrmzZunaDQa5fr161myvh07digajUbZv39/hnEpKSlK165dlSpVqihHjx5NtXzp0qWKRqNRpk2bliXlMnYajUYZNWqUOn3s2DFFo9EomzdvfiPbv379eqoy5ISGDRsqXbt2zdEyZMTQ8sXFxSmenp5KrVq1lHPnzqVafv36dcXNzU1p1KiRkpCQkB1FfSvI7cG3xLp164iIiGDMmDFUrFgx1fJRo0ZRqFAh1q9fnwOly166wUl1L49Mz759+zhx4gS9e/emTp06qZb36tWLcuXKsW3bNp49e5YtZRUip6xYsYLr168zcuRIHB0dUy0vVaoUgwYN4tq1axw7diwHSvhmSNJ6S+zcuZP33nuPFi1apLk8X758bNiwgW3btunNP3nyJH5+fmi1WrRaLd26dUs1YrOHhwe+vr6p1vnyfA8PDyZMmMD27dtp0aIFVapUoXHjxqxdu1aNGT16NAsWLADA09MzzfW+6O+//+azzz7DxcWFqlWr4uPjo77qG8DX11d9Q2q3bt1wcHBId107d+4E0LvV9LLvvvuOQ4cOkS9fPnXeTz/9RNeuXalRowZOTk54eHgQGBio1/6VkJDAl19+iaenJ05OTtSvX5/JkycTHR2tt/7bt28zcuRIPvroI6pUqUKbNm3YsWOHXoyiKCxYsIAmTZpQpUoV6tSpw4gRI7h161aG++rlNi3d9MWLFxk+fDg1a9ZEq9XSv39/bty4ke56bty4oe7HrVu3pmoni42NZfLkydSuXZtq1arRvXt3/v77b711pKSkEBQURNOmTXFycsLNzY2pU6dm2Vtqde1rW7duxdvbmypVqqi/A0O3fe/ePb744gv1O6tRowbdunXj999/V2McHByIjIzkxIkTanuebv/s2LGDmTNnUqdOHbRaLZ999hkPHjzg7NmzdOrUCWdnZ5o0acKuXbsyvW909Tt69Ki6r52dnenevTsXL17MsHzp+fHHHylevHiGt3abNm3K2rVr1ZM6XV1XrlxJ586dcXJy0nuh58aNG2ndujVVqlTho48+Yvjw4Xq/Ld3nXx5z8+X5uunNmzczY8YMXF1dqVmzJoMGDcry9kQZ5f0toCgKFy5coHr16hkOPlm2bFm96f379zNgwAB1tHV4/iP08/Nj3rx5eHp6ZrosR44cUQ/yxYoV44cffuCLL76gVKlS1K9fn44dOxITE8O+ffsYM2YMH374YbrrOnv2LN26dcPKyooePXpgaWnJ9u3b6d+/PxMmTKBLly707dsXe3t7fvjhB/r27Uu5cuXSXd9ff/2Fra1tqvdbvejlNraNGzcyfvx4PDw8CAgIIDExkX379rF8+XLee+89BgwYAMAXX3xBSEgI3bp1o3Tp0vz777+sXbuWq1evEhQUBDx/yWeHDh1QFAVfX18KFSrE/v37GTFiBHfv3qVnz54ALF68mG+//ZYuXbrg4ODAjRs3WLVqFefOnSMkJAQzMzODvw+Afv36Ub58eYYOHcr169f5/vvvuXPnDps2bUoz3tramsDAQEaOHKm2J5UvX159M+9XX32Fg4MDAwcO5M6dO6xYsQJ/f39CQ0PVZD9u3Di2bdtG27Zt8fPz4/LlywQHB3Pq1CmCg4P1Bo99HV988QWtW7emQ4cO6mjhhmz72bNndOnShSdPntClSxdsbGy4cuUKwcHB9O7dm8OHD2NlZUVgYCDTp0+nSJEi9O3bV++9UV999RXFixdnwIABXLp0ibVr1/Lw4UPCw8P5+OOP8fb2ZtWqVeqVTZkyZTK9b8aPH8/777/PZ599RnR0NMuWLaNXr14cPHgQc3PzDMv3oqioKC5fvkzLli0zbNO2srJSX6Hyom+++Yb69evj7e2tlm/mzJkEBQVRu3ZtRo4cyd27d1mzZg2//PILGzdu/E/t1QsWLCAlJQV/f3/i4uJYuXIlp0+fZseOHalepPqf5fDtSaEoSlRUlKLRaJShQ4ca/JnExETF3d1dqV+/vvLkyRN1fnR0tOLm5qa4ubmp97XTu2f+8vyGDRsqDg4Oem1qd+/eVRwcHJRhw4ap8wxt0+rQoYNSrVo15datW+q8Z8+eKW3btlWqVq2qREVFKYqiKJs3b1Y0Go1y7NixDNfn7Oys+Pj4ZBjzsqZNmyodO3ZUUlJS1Hm6fdeyZUt1XtWqVZXJkyfrfXbOnDnKxx9/rMTExCiKoiijRo1SatWqpdy5c0cvbtiwYYqTk5Ny//59RVEUpVmzZkrv3r31YoKDg5VWrVopV69eTbesL+8H3fSAAQP04iZMmKBoNBolIiIiw7qn16bVtm1bJTExUZ0/f/58RaPRKL/88oteXHBwsN76jhw5omg0GmXlypUZblcnozYt3TZe/l0auu2dO3cqGo1G+fnnn/XigoODFY1Go+zZs0ed9/LvXFcud3d3JS4uTp3frl07RaPRKOvWrVPnHT16VNFoNMoPP/yQqfLp4tq1a6ckJSWpcUuWLFE0Go0SFhaWbvnScu7cOUWj0SgzZ85Mtezx48dKVFSU3j/db1ZXVy8vL712rkuXLikODg5K//799f42/vjjD8XBwUEZPHiw3ufnzZunt82X5+umX/571+2/r776KsP6ZYbcHnwLmJo+/xp0b0g1xPnz57l9+zZdunTRe31AwYIF6dq1K3fu3OHcuXOZLou9vb1em1rx4sUpVqwY9+/fz9R67t+/z5kzZ2jdurX62hZ4/uJIf39/nj17xi+//JKpdZqammZqHwHs2LGDpUuX6p2dRkVFUbBgQb0u+yVKlGDXrl1s2bKFx48fAzBkyBA2b96MpaUlKSkphIaG4uLigrm5OQ8ePFD/NW7cmISEBI4ePaqu6/jx43z//ffqfuvUqRPbt2/Hzs4uU+UHaNasmd50pUqVADL9neg0bdoUc/P/3WTR9dTUrW/v3r2YmJhQv359vXpWrlyZ4sWLq+/Uygr16tXTmzZ0282bN+fXX3/V+/yLt3sNeRzDzc1N7zay7k6Gl5eXOk93tXH37t1MlU+ncePGelfWuu/u3r17ryzfi3SPuyhpdPbu3r07tWvX1vs3ZcoUvZiPPvpI7y7OgQMHUBSF3r176/1tODs7U7duXQ4dOkRSUlKmygjQqlUrvb/3OnXqoNFo1BdpZgW5PfgWKFSoEHny5OHBgwcGf0Z33zmtV13rbrHdvHkz011fra2tU82zsLDI9DNikZGR6ZZP93r6zD5TUrx48UztI4A8efLw22+/ERISQnh4ONeuXSMqKgoAW1tbNW7SpEkMGTKEMWPG8Pnnn1OtWjW8vLxo164dBQoU4OHDhzx58oTQ0FC9NrkX6dqsRo4cSb9+/Zg2bRrTp0/H0dERDw8PfHx8KF68eKbKD6S6raLrzp/ZBK7z8nesO3DrOsRcu3YNRVFo0KBBmp9/VYeZ1ylLZrZtYmLC0qVLOX36NNeuXePatWtqHQz5vRYtWlRvWpfIX5yvO6HUJYvM7puX66f77jL796S7JZ7WicqkSZP02tN69OiRKublcmR0/ChfvjxhYWE8fPgwU2UEqFChQqp5ZcuW5ciRI5leV3okab0FTExM0Gq1nDt3jqSkJL2z4BfNmTOH69evM2bMmDTPuHR0y171cra0Dnq6P9LXlVH5dH+wmX15nFarZevWrdy9ezfddq21a9fyyy+/MGzYMMqXL8/XX3/N0qVLqVy5MtWqVaN169ZotVqmTJmi1zGidu3aHDx4UP139OhRpk+fzsqVK9myZYu6r5o0aUKnTp3S3Hbp0qUBqFixInv27OHIkSMcPHiQI0eOMG/ePFauXMn69evVpG2orPpODF1fSkoKlpaWaoebl2VVexaQqn3P0G1HRkbSsWNHYmNjqVevHs2bN6dSpUooikL//v0N2nZ6f2cZtRlldt9k1XdnY2ODra0tJ06cSLWsatWqr/z8y/vZ0L/P+Pj4DGNeltbfdHJycpb+hiVpvSW8vLw4ceIEu3btolWrVqmWP3v2jE2bNpGcnEzhwoXVq4Tw8PBUsREREQDqZbqpqWmqkSKSkpJ4+PDhf7pdZYjMlM9QXl5ebN26lc2bN6sdT16UkpLCxo0buXTpEpMnTyYyMpKlS5fSunVrAgMD9WJfPGNNSEjgwoULlChRghYtWtCiRQtSUlJYsWIFgYGB7Ny5k86dO5M/f36SkpJSdbe/efMm58+fJ3/+/CQnJ3Px4kWsrKzw9PRUO8Ps2rWLoUOHsnHjRkaPHp2per9ptra2hIWF4eTkRMGCBfWW7dmzh8KFC+f4thcsWEBUVBS7d+/W66D0448/ZlvZMlO+7ODt7c3ixYvZvXt3qlvGmaW77RkeHo6zs7PesoiICN577z0KFSpEXFwcQKrjR3q3N69du5Zq3tWrV1N1Insd0qb1lujYsSO2trbMnDmTf/75R29ZcnIykyZN4v79+/Tq1Ys8efLg6OhI8eLFCQ4O1rs1EBMTw7p16yhevDhOTk4AFCtWjIiICL1nlw4cOJDuWdSrvHzLJC267e/YsYPbt2+r8xMSElixYgUWFhbUrVs3U9v19PSkWrVqLF26NM2hjubPn8+FCxfw8fGhWLFianf1l29ZHD58mCtXrqj37B89ekTHjh1ZsmSJXh11bT2mpqaYm5vj7u7O4cOH9bosA8yYMYP+/fvz8OFDkpOT6datG9OmTdOL0R0YsvqqKSOmpqb/aegv3TBLixYt0pt/4MABBg0alK2JwdBtP3r0iPz586s9DuH5b0v3HOOLdxH+6354nfJlhqHl69OnD3Z2dkyaNImTJ0+mWh4fH8/XX39t0DYbNmwIPH9E5MW/47/++otffvmF+vXrY2JiQuHChTE3N+fChQt6n9+9e3ea692+fbve8ejw4cNcunSJxo0bG1QuQ8iV1lsib968LFiwgE8//ZT27durz648evSIn376iQsXLtC0aVP1fnWePHn4/PPPGTJkCO3ataN9+/YAbNq0ibt37zJv3jz1ANmyZUumTJlCz549adWqFVevXmXDhg16bTqZobs/vmzZMtzd3dPtWj9+/Hi6d+9O+/bt6dy5M5aWluzYsYO//vqL8ePHpzpTNcTs2bPx8/OjR48eNGnSBBcXF+Li4jh06BC//fYbLi4uBAQEAM+TVcmSJVm8eDHx8fGUKFGCs2fPsnXrVvLmzcvTp0+B5+0F3t7erFu3jri4OLRaLY8ePWLNmjUUK1ZMPasNCAjg+PHjdOnShS5dulCyZEkOHTrEwYMH6dixo9r939fXl0WLFtG/f3/c3Nx49uwZP/zwA/nz56ddu3aZrvN/ZW1tzYkTJ9iwYUOqDg8ZqV+/Pp6engQFBXHjxg3q1KlDZGQka9eupWTJkvj7+2dbmQ3dtru7OwcOHKBPnz40bdqUJ0+esG3bNvVMX/fdwvP9cPHiRdatW0etWrX0Ol9kV/ky4+XypdUuBPDee+/x3XffMXDgQLp164a7uzuurq5YWlpy+fJldu7cyb1793B0dEyzXetFH374Ib6+vqxevZoePXrQqFEj7t27x+rVqylYsCDDhw8HIH/+/Hh6erJnzx7Gjh2LVqvl2LFjnD59Os1bgU+ePMHHxwcfHx8ePHjA999/T/ny5enevXum90t6JGm9RSpXrsz27dtZuXIlP//8M7t27UJRFBwcHJg2bRoff/yx3v32Jk2aEBQUxMKFC/n2228xNzfH2dmZL7/8Uu9ZjU8++YRHjx6xadMmpkyZQsWKFVmwYAFBQUH/adDbFi1asHfvXrZs2cKJEyfSTVparZbg4GDmzZtHUFAQKSkpVKxYkW+//fY/jz1na2vLxo0bWbdundpulJCQQLly5Rg9ejRdu3ZV/5gsLCxYunQpM2bMYNWqVSiKgp2dHWPHjiUpKYkvv/ySc+fO4eTkxJQpUyhdujQ7d+5k586d5M+fn9q1azN06FA1SdvZ2bFhwwbmzZvHhg0biI2NpXTp0owZM0bvIetBgwZRuHBhNm/ezMyZMzEzM6N69erMmjUr0+1ZryMgIICvv/6aKVOmMGXKFD744AODPmdiYsI333zDsmXL2LZtGwcPHsTa2prGjRszePDgLBkX8XW33alTJx4/fszGjRuZOnUqxYoVo1q1aixYsIBOnTpx7Ngx9SHagQMHMnHiRKZNm0b//v3x9vbO9vJlxsvlSy9pwfNODRs3bmT79u2EhISwbNkyoqOjKVq0KDVq1KBVq1Z4eHgYND7puHHjsLe3Z/369cyYMYNChQrh5eXFoEGD9E5oJ0+ezHvvvUdoaCi7d++mdu3arF69mpYtW6ZaZ48ePYiJiWH+/PlYWFjQqlUrhg0bRv78+TO9X9IjA+YKIYR4LTdu3MDT05MBAwYwcODAbN2WtGkJIYQwGpK0hBBCGA1JWkIIIYyGtGkJIYQwGnKlJYQQwmhI0hJCCGE0JGkJIYQwGpK0hBBCGA1JWkIIIYyGJC0hhBBG4/8Ah28aBtDteS8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(simulations_fast, bins=30)\n",
    "#plt.axvline(0.57, color='red', linestyle='dashed', linewidth=2)\n",
    "plt.title('Approximate Sampling Distribution')\n",
    "plt.ylabel('# of Simulations')\n",
    "plt.xlabel('Count of Cases in the Treatment Group')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our simulation, we repeated the process of random assignment to the treatment group 500,000 times.\n",
    "Indeed, we found that on average there were 234 cases of disease in the treatment group.\n",
    "The histogram of outcomes shows that not one of the 500,000 simulations had as few as 117 cases or fewer.\n",
    "It would be an extremely rare event to see so few cases of COVID-19, if in fact the vaccine was not effective. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.754385504532372"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.std(simulations_fast)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vaccine Efficacy\n",
    "\n",
    "Vaccine Efficacy (VE) is measured by comparing the risk of disease among vaccinated and unvaccinated persons: \n",
    "\n",
    "$$ \\frac{\\textrm{Risk among unvaccinated group} − \\textrm{Risk among vaccinated group}}\n",
    " {\\textrm{Risk among unvaccinated group}},$$ \n",
    " \n",
    "where, e.g., the risk among unvaccinated is the proportion of unvaccinated who contracted COVID. \n",
    "\n",
    "Since the two groups had the same number in each, we can compute the efficacy as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6666666666666666"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(351 - 117) / 351"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Centers for Disease Control sets a standard for VE of 50% when deciding whether to adopt a new vaccine. This would be equivalent to how many cases in the treatment group?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "156.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "468/3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see from the histogram that none of the simulations yielded as few or fewer than 156 cases in the treatment group. Furthermore, when it came to preventing severe cases of Covid, the J&J vaccine was over 80% effective, and no deaths were observed in the treatment group. \n",
    "\n",
    "The Johnson & Johnson vaccine received emergency approval in spring 2021, and expects to receive full approval by the FDA in fall 2021. \n",
    "The Mayor of Detroit retracted his original statement, saying \n",
    "“I have full confidence that the Johnson & Johnson vaccine is both safe and effective.” "
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
